Slip flow over structured surfaces with entrapped microbubbles 
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On hydrophobic surfaces, roughness may lead to a transition to a superhydrophobic state, where 
gas bubbles at the surface can have a strong impact on a detected slip. We present two-phase lattice 
Boltzmann simulations of a Couette flow over structured surfaces with attached gas bubbles. Even 
though the bubbles add slippery surfaces to the channel, they can cause negative slip to appear 
due to the increased roughness. The simulation method used allows the bubbles to deform due 
to viscous stresses. We find a decrease of the detected slip with increasing shear rate which is in 
contrast to some recent experimental results implicating that bubble deformation cannot account 
for these experiments. Possible applications of bubble surfaces in microfluidic devices are discussed. 

PACS numbers: 83.50.Rp, 47.55.D-, 47.11.-j 
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The no-slip boundary condition states that the fluid 
velocity at a fluid-solid interface equals to that of the 
solid surface. Although this boundary condition has been 
proven valid at macroscopic scales, it has no microscopic 
justification and is not fulfilled generally 1]. Its valid- 
ity was debated already in the early days of fluid me- 
chanics and due to recent developments in microfluidics 
the interest in violation of the no-slip boundary condi- 
tion has revived [2j. In microfluidics, several experiments 
have found fluid slip at the boundaries of the flow chan- 
nels 0, 0j- As the slip length (defined below) has typ- 
ically a magnitude measured in nano- or micrometers, 
the appearance of slip does not have noticeable ramifi- 
cations in macroscopic flows. However, in microfluidic 
devices with large surface-to- volume ratio, surface prop- 
erties may dramatically affect the flow resistance. The 
possibility to engineer slip properties in a controlled way 
is therefore crucial for microfluidic applications. 

The slip at fluid-solid boundaries can be quantified by 
Navier's boundary condition, which states that the slip 
velocity is proportional to the velocity gradient, i.e., v s — 
b(dv/dz) at the surface z = z - Here, b is the slip length 
which is the distance below the surface where the velocity 
linearly extrapolates to zero. Higher slip length means 
larger slip and lower friction at the boundary. 

Mechanisms behind boundary slip include surface 
roughness and structural characteristics of roughness [!, 
6], roughness- induced dewetting on hydrophobic sur- 
faces ti_|8|i__ dissolved gas and bubbles on the sur- 
face [3, [lfj, [ll( , as well as wetting properties [12j . Usu- 



ally, roughness decreases the slip length due to increased 
dissipation and the roughness induced slip is just an ar- 
tifact [f|. Theoretically, it was shown that even slippery 
surfaces, if rough enough, can provide no-slip bound- 
aries [l3j. However, if surfaces are hydrophobic, rough- 
ness may increase the slip due to a transition to a su- 
perhydrophobic (fakir or Cassie) state Here, liquid 
cannot enter between roughness elements but stays at 
the top of them. Thus, gas bubbles or layers are formed 



thereby lubricating the flow due to a reduced liquid-solid 
contact area. By using a surface patterned with a square 
array of cylindrical holes, Steinberger et al. found that 
gas bubbles may also cause an opposite effect, i.e., slip 
is reduced if microbubbles are present in the holes [1QJ. 
Numerically, they found even negative slip lengths for 
flow over such a bubble mattress. Negative slip means 
that the effective no-slip plane is inside the channel, i.e., 
the bubbles increase the flow resistance. These mecha- 
nisms are related to so-called effective slip and should be 
distinguished from the (smaller) intrinsic slip on smooth 
surfaces. Another peculiarity in the boundary slip is the 
shear-rate dependence observed in some experiments but 
not in the others Q. The question if the shear-rate de- 
pendence is a true property of slip is still to be answered. 

As the different mechanisms behind the slip phe- 
nomenon are strongly intertwined, the experimental 
study of a single mechanism is a complicated task. There- 
fore, numerical simulations are attractive as they provide 
a controllable way to change a single property of fluid 
or surface while keeping the others unchanged. Most 
computer simulations so far have been performed us- 
ing molecular dynamics 0, EE EE [l7j ]. For computa- 
tional reasons molecular dynamics is limited to length 
scales of tens of nanometers and time scales of nanosec- 
onds, which do not comply with the experimentally 
relevant scales. Therefore, mesoscopic lattice Boltz- 
mann (LB) simulations have recently been applied to 
study flow in microchannels or along hydrophobic sur- 
faces [E H, [H, EE Ell- This method allows to reach 
experimentally relevant scales and preserves those inter- 
actions needed to describe the underlying physics. 

Our simulations utilize the multiphase LB model by 
Shan and Chen ,20] . Dynamics of the method is governed 
by a discretized Boltzmann equation 

/,(r + c, t + 1) - fi(r, t) = ~ [Mr, t) - Jf (r, f)] , (1) 

where /» is a distribution function describing the proba- 
bility to find a particle at position r at time step t, mov- 
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FIG. 1: (Color online) Slip length & as a function of protrusion 
angle <p. A unit cell of each array is shown in insets and 
corresponding results are given by triangles (rhombic array), 
diamonds (rectangular array), and circles (square array). The 
inset in the top-left corner shows the definition of ip. 



ing in lattice direction Cj. We use a three-dimensional 
lattice with 19 discrete velocities. The right-hand side 
of Eq. Q] models the relaxation of the fi towards a lo- 
cal equilibrium due to collisions among the particles on a 
time scale given by the relaxation time r. The mean-field 
interactions between particles are given by a force 



F = GMA(r) J>V(r - 



(2) 



where ip = 1 — exp (— p/po) is an effective mass (p is 
the fluid density, po a reference density), Gb tunes the 
strength of the interaction, and ij's are weight factors 
for different lattice directions. This force term leads 
to a non-ideal equation of state with pressure P — 
c 2 s p + TjC^Gbip 2 , where c s is the speed of sound, and it 
enables simulations of liquid-vapor systems with surface 
tension. To model the wetting behavior at fluid-solid sur- 
faces, a similar interaction is added between the fluid and 
solid phases, and the contact angle is tuned by setting a 
density value p w at the boundaries [21 1. 

In this article we investigate liquid slip in Couette flow, 
where the flow is confined between two parallel walls. 
One of the walls is patterned with holes and vapor bub- 
bles are trapped to these holes. Steinberger et al. [lo| 
presented finite-element simulations of flow over rigid 
"bubbles" by applying slip boundaries at static bubble 
surfaces. The LB method allows the bubbles to deform 
if the viscous forces are high enough compared to the 
surface tension. We are also interested in how surface 
patterning affects the slip properties of these surfaces, 
and how bubbles could be utilized to develop surfaces 
with special properties for microfluidic applications. 

In our simulations, the lower wall is static and has 
the topographical patterning whereas the upper one is 



smooth and moved with velocity uq. The distance be- 
tween walls is d = 1 /xm (40 nodes) in all simulations, and 
the area fraction of holes 0.43 unless stated otherwise. 
The system boundaries are periodic and a unit cell of the 
regular array is included in a simulation. To trap bub- 
bles to holes, some heterogeneity is needed at the edges of 
holes in order to pin the contact line. To this end, we use 
different wettabilities for boundaries in contact with the 
main channel and with the hole. The protrusion angle ip 
(see Fig. [T] for definition) is varied by changing the bulk 
pressure of the liquid. A similar technique to form bub- 
bles on structured surfaces was used experimentally by 
Bremond et al. 22]. The effective slip length can be cal- 
culated from the shear stress a — pdv/dz acting on the 
upper wall, which is obtained from the no-slip boundary 
condition imposed at the fluid-solid boundaries. Thus the 
effective slip length reads as b — puo/cr — d, where p is the 
dynamic viscosity of the liquid, b is measured from the 
top of the structured surface and Gb is chosen such that 
the density ratio between liquid and gas is 22. This ratio 
is too small for a realistic description of gas bubbles in a 
liquid. Also, the interface between both phases is of finite 
width causing the resistance in the vapor phase. These 
limitations of multiphase LB models do not influence the 
qualitative insight obtained from our simulations. 

In order to study the effect of a modified protrusion 
angle and different surface patternings, we use three dif- 
ferent arrays of bubbles, i.e. , a square array, a rectangular 
array where the distance of bubbles is larger in one di- 
rection than in the other, and a rhombic array. These 
surfaces have cylindrical holes with radius a = 500 nm 
and the area fraction of the holes is equal in all cases. The 
shear rate is such that the Capillary number Ca — 0.16. 
The Capillary number is the ratio of viscous and sur- 
face forces, i.e., Ca = paG s /j, where G s and 7 are the 
shear rate and surface tension, respectively. Snapshots 
of simulations are shown in Fig. and the slip lengths 
obtained are shown in Fig. The observed behavior is 
similar to that reported in [lfj, where a square array of 
holes was studied. In particular, we observe that when 
the protrusion angle is large enough the slip length be- 
comes negative. We also see that the maximum of the 
slip length is obtained when the protrusion angle equals 
zero. Since the area fraction of the bubbles is the same in 
all three cases, the results clearly indicate that slip prop- 
erties of the surface can be tailored not only by changing 
the protrusion angle but also by the array geometry. 

Next, the shear-rate dependence of the slip length is 
investigated. As the shear rate and thus the viscous 
stresses grow the bubbles are deformed (see Fig. and 
the flow field is modified. Results are shown in Fig. [3^, 
where the slip length is given as a function of the Cap- 
illary number for three different protrusion angles. We 
also calculate the Taylor deformation D = (i — a)/(i + a) 
of the bubbles by fitting an ellipse to the bubble inter- 
face. Here, a and I are the minor and major axes of 



FIG. 2: (Color online) Snapshots of simulations of bubbles on structured surfaces. Shown are a square array of bubbles with 
Capillary number (a) Ca — 0.02 and (b) Ca — 0.4, and (c) a slot with a cylindrical bubble. In each case a unit cell is shown. 
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FIG. 3: (Color online) (a) Slip length as a function of capil- 
lary number for a square array of bubbles with three different 
protrusion angles, tp — 63°, 68°, and 71° (from uppermost to 
lowermost). In the inset shown are cross sections of liquid- 
gas interfaces for four capillary numbers, (b) Taylor defor- 
mation as a function of Capillary number for the bubble with 
tp — 71°. These Ca values correspond to shear rates of the 
order of 10" 6 - 10~ 7 s _1 . 



the ellipse, respectively. A slightly superlinear relation 
between D and Ca is observed (Fig. [3b) in accordance 
with a two-dimensional case in Ref. [23j| . We find that 
increasing shear rates cause the slip length to decrease. 
This behavior is contradictory to those found in some 
experiments using surface force apparatuses (see. e.g., 
Ref. [HI), where a strong increase in the slip is observed 
after some critical shear rate. This shear-rate dependence 
has been explained, e.g., with formation and growth of 
bubbles 0, . However, one should notice that these 
experiments are dynamic in nature while we simulate a 
steady case. In our simulations there is no formation 
or growth of the bubbles but we determine the slip for 
given bubbles of given size. Our results indicate that 
the deformation of the bubbles and the changes in the 
flow field thereby occurring cannot be an explanation for 
the observed shear-rate dependence. On the other hand, 
our results are consistent with [H , where it is shown that 
smaller roughness leads to smaller values of a detected 
slip. In the present case, the shear reduces the average 
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FIG. 4: Slip length as a function of hole area fraction. Filled 
circles denote values for the flow parallel to the slots and 
triangles relate to the perpendicular direction. Correspond- 
ing open symbols are for a surface with rectangular holes 
(tp = 72°). The inset shows the normalized slip length as 
a function of the solid-area fraction, where the dashed line is 
the theoretical prediction for small solid-area fractions. 



height of the bubbles and thus the average scale of the 
roughness decreases as well. 

As seen above, the slip properties of a bubble mattress 
can be tailored by changing the surface patterning. Next, 
we consider a surface which has slots with a width of 
1 /xm, and cylindrical bubbles protruding to the channel 
with an angle tp = 72° (Fig. . We apply shear parallel 
and perpendicular to the slots, while the area fraction of 
the bubbles (0) is varied by changing the distance be- 
tween the grooves (L). According to our results (Fig. [4]), 
b strongly depends on the flow direction and even the 
qualitative behavior changes. When the flow is parallel 
to the slot a positive slip is observed, but for perpen- 
dicular flow the slip becomes negative. This observation 
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can be explained by means of theoretical predictions of 
Richardson, who showed that slippery surfaces lead to 
(macroscopic) no-slip boundaries if rough enough [l3j |. 
Our results support these predictions. In the case with 
flow parallel to the bubbles the streamlines are straight 
and the flow does not "see" any roughness. However, in 
the perpendicular direction roughness is caused by the 
bubbles and negative slip is observed. The inset of Fig. [4] 
depicts the scaling behavior of the slip length as a func- 
tion of the solid-area fraction for flow perpendicular to 
the slots. We see that for small solid-area fraction the 
slip length obeys the scaling law b ~ — Llog(l — 0) as 
was recently predicted in Ref. [f|. Anisotropic behavior 
has been observed in the case of flat surfaces composed 
of stripes with no-slip and perfect-slip boundary condi- 
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261 ]. Our simulations differ from these studies 
as our bubbles protrude to the channel thus leading to 
a larger effect than for a flat surface. In order to under- 
stand less idealized surface patternings, we study flow 
over rectangular holes. The aspect ratio of the holes is 
chosen such that the longer side is three times the smaller 
one. We observe similar qualitative behavior (Fig. 2]), but 
the difference in slip length between the two flow direc- 
tions is less pronounced. Obviously, by changing the as- 
pect ratio of the hole, the anisotropic behavior of the slip 
can be tailored. Due to the striking difference between 
the slip properties in the two perpendicular directions, 
we believe that this kind of surfaces may find applica- 
tions in microfluidic devices. Anisotropic surfaces could 
be used to suppress the flow in an unwanted direction 
while enhancing it in another one. The suppressing be- 
havior is further amplified by the shear-rate dependence, 
since the negative slip is growing with the shear rate. 

To conclude, we simulated Couette flow in a mi- 
crochannel where one of the walls is patterned and mi- 
crobubbles are attached to the pattern. We found that 
the slip properties of the surface can be tailored by chang- 
ing the hole array and that such a surface with bub- 
bles may yield negative slip, i.e., increased resistance to 
flow, if bubbles are strongly protruding to the channel. 
Our results can be qualitatively compared to previous re- 
sults [10| | , but overcome their limitation of a static liquid- 
vapor interface. This allowed to study the influence of 
the shear rate on the deformation of the interface and 
it's effect on the measured slip. We found that the slip 
decreases with increasing shear rate demonstrating that 
shear induced bubble deformation cannot explain recent 
experimental findings where slip increases with increasing 
shear rate [24| . We proposed a special surface patterning 
which can be used to produce surfaces where the slip is 
positive in one direction and negative in the perpendic- 
ular one. Such a surface might be useful to construct 
microfluidic devices with tunable flow throughput which 
could be controlled by adding bubble ridges parallel or 
perpendicular to the flow. In addition, we have shown 
that the throughput could be tailored by tuning the bulk 



pressure, i.e., the protrusion angle, or the shear rate. 
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